Monte Carlo simulations on the Lefschetz thimble: taming the sign problem 
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We present the first practical Monte Carlo calculations of the recently proposed Lefschetz thimble 
formulation of quantum field theories. Our results provide strong evidence that the numerical sign 
problem that afflicts Monte Carlo calculations of models with complex actions can be softened 
significantly by changing the domain of integration to the Lefschetz thimble or approximations 
thereof. We study the interacting complex scalar field theory (relativistic Bose gas) in lattices of 
size up to 8 4 using a computationally inexpensive approximation of the Lefschetz thimble. Our 
results are in excellent agreement with known results. We show that — at least in the case of the 
relativistic Bose gas — the thimble can be systematically approached and the remaining residual 
phase leads to a much more tractable sign problem (if at all) than the original formulation. This is 
especially encouraging in view of the wide applicability — in principle — of our method to quantum 
field theories with a sign problem. We believe that this opens up new possibilities for accurate Monte 
Carlo calculations in strongly interacting systems of sizes much larger that previously possible. 



Introduction — Many important physical systems are 
characterized by complex actions, when formulated in 
terms of a path integral. But, if the action S is not 
real, then e~ s is not positive semi-definite and it can- 
not be interpreted as a probability distribution. In these 
cases, Monte Carlo calculations are not applicable di- 
rectly. This is the so called sign problem. Many tech- 
niques have been proposed to overcome this problem, 
with important partial successes, but the sign problem 
is still unsolved for a variety of important physical sys- 
tems and parameter values, such as lattice QCD at high 
baryonic density [T], or with a 9- vacuum [2], real-time 
quantum field theories [3] , the electron structure calcula- 
tions 0HB] , the repulsive Hubbard model [7] , the nuclear 
shell model [S] or polymer theory [S], to mention only 
some of the most famous problems. In this context, any 
new idea that could improve our chances to simulate any 
of these models on larger lattices than are feasible today 
would be extremely valuable. 

In a previous work [101 1 1 1 j , we argued that it may be 
possible to control the sign problem by reformulating the 
associated quantum field theory on a Lefschetz thimble. 
The Lefschetz thimble, associated with a saddle point <fi, 
is defined as the hypersurface formed by the union of 
all paths of steepest descent (SD) of the complex action 
ending in that saddle point <p. Both the Lefschetz thim- 
ble and the saddle point are constructed in an enlarged 
space obtained by complexifying each field component. 
We showed that, in many cases of interest, this reformu- 
lation has the same symmetries and perturbation theory 
as the original theory [TU]. Thereafter, appealing to uni- 
versality we argued that the reformulation has the same 
physical content as the original theory. 

The benefit of this reformulation is that the action 
on the Lefschetz thimble has a constant imaginary part, 
which can be set to zero without any loss of generality. 



Thus e~^ s ^ can now be interpreted as a probability dis- 
tribution in Monte Carlo sampling. Since, the Lefschetz 
thimble defines a curved integration domain, there can, 
in principle, be an additional residual phase coming from 
the Jacobian of the transformation. However, we will ar- 
gue later that this residual phase, if at all present, will 
result in a very mild growth of stochastic noise. 

In this work, we apply our method to the interacting 
complex scalar field theory describing a relativistic Bose 
gas at finite chemical potential. This model is one of the 
simplest non-trivial examples whose sign problem shares 
many features with the more complex systems mentioned 
above. Also, in common with lattice QCD, it displays the 
Silver Blaze phenomenon |12j . i.e., the independence of 
the physics on the chemical potential up to some (finite) 
critical value. This feature is not accessible to standard 
Monte Carlo treatments due to the sign problem. Quite 
importantly, this model has been solved through alterna- 
tive methods [T5lU6| . and as such provides the ideal test 
bed for new methods, like ours, for studying the physics 
of strongly interacting systems. 

In the context of Monte Carlo methods, modifications 
of the domain of integration had been proposed already in 
[T71 [TS] . But those deformations were limited to shifts 
of the contour in the imaginary direction. For many rele- 
vant theories, including those considered in [10], the shift 
is zero, and more general transformations are necessary, 
to reduce the sign problem. Morse theory jTQl [THl HO] 
identifies the Lefschetz thimbles as the appropriate con- 
tours of integration in the more general cases. 

Formulation of the model on a Lefschetz thimble - 
The model is denned by the following continuum action: 

S = J d^xim 2 + (m 2 - M 2 )|0| 2 + W - + A|0| 4 ], (1) 
where <fi(x) is a complex scalar field, j v := <jf <9„0 — </>9„0* 
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and n is the chemical potential. In this model (as in 
QCD) the density (n) = ydlnZ/dfj, is expected to be 
zero up to a critical point. But, this phase transition is 
hidden in the standard Monte Carlo method because of 
the strong sign problem which appears as soon as /i ^ 0. 

To formulate and simulate the relativistic Bose gas on 
a Lefschetz thimble [TO], we need to discretize the sys- 
tem defined by Eq. ([!]) and extend the action S holo- 
morphically. This is done by complexifying both the 
real and imaginary part of the original complex fields 
4>x = ^{4>l,x + i<t>2,x), a S 4>a,x = 4%$ + i9a}x, a = 1,2, 
which leads to the action in d dimensions [2"T] : 



help of Langevin dynamics using an algorithm described 
in [TQl [11], that we review here. First, let us assume to 
know a starting configuration (f> G j7o, together with a 
set of configurations 0(/cAr) € j7o, with k = 1, . . . , N T , 
that represent the path of SD connecting tfi = 0(0) with 
the configuration 0(r — AtN t ). Let us assume that 
0(r) is sufficiently close to g i o b, so that the action S 
can be approximated by its quadratic expansion around 
'/'glob- Second, we generate a Gaussian noise rjj, where 
j = 1 . . . 2N is a multi-index that stands for (R/I, a, x), 
and we project it as follows: 
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(e is the 2 dimensional anti-symmetric Levi-Civita sym- 
bol). The observables are defined as: 
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where the integration domain j7o is the Lefschetz thimble 
[THEO] attached to g i o b- The configuration g i o b is the 
global minimum of the real part of the action Sr = 3i{S}, 
when restricted to the original domain M. 2V . More pre- 
cisely, j7o is the manifold of real dimension N = 2V, 
defined as union of all the curves of SD for Sr, i.e., the 
curves that are solutions of: 
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(4) 



and that end in g i o b for r — > oo. 

In presence of spontaneous symmetry breaking (SSB), 
the global minimum g i o b is degenerate. But, the whole 
procedure can be defined by introducing an explicit term 
of symmetry breaking: h J2 X a 4>x,a-, where h is a real con- 
stant, that selects a specific minimum [TJJ (that can be 
computed also analytically). Since h is real, the global 
minimum g i o b of Sr is also a stationary point of the 
imaginary part of the action Sj, and hence the thimble 
is well defined. Physical results are obtained by extrapo- 
lating to h —> 0. 

Aurora Monte Carlo algorithm for sampling the thim- 
ble — It is possible to generate field configurations on the 
Lefschetz thimble with weights given by e _Sf? with the 
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Pj,kVk, 



(5) 



where the 2N x 2N matrix P of rank N is defined in 
terms of the Hessian matrix H as: 



P 



H 



1. 



and 



H = d 2 S R [4 gloh }. (6) 



Then we normalize the noise vector as: 



rj = r- 



\V 



(7) 



where r is a random number distributed according to the 
iV-dimensional \ distribution. This produces a Gaussian 
noise on the tangent space to j7o computed in </> g i b- We 
call such linear vector space Q$. Third, we transport the 
noise from s = r along the path of steepest ascent (SA) 
to s = by integrating the ordinary differential equation 
(ODE): 



ds 



v'(s)=yn' k (s)A(s) K 



k 



J' 



(8) 



where A(s)k t j is the Iwasawa projection of dkdj S r[4>(s)] 
|11) . This ensures that the noise remains tangent to J§. 
Fourth, we use the evolved noise r?'(0) to generate a new 
configuration as: 



^ = 4>j - At 



6 S R [<I>] 



In the limit At — > this simulates Langevin dynamics on 
the thimble. For At > 0, <fi'(0) will move away from the 
thimble of order (At) 2 . To correct this, the fifth step con- 
sists in following the path of SD from (j>'(0) for a length r 
leading to the configuration 0'(t). Assuming that the ac- 
tion at 4>' ( r ) can be approximated with its quadratic part 
(otherwise, we extend r), we ensure that 0'(t) belongs to 
the thimble by projecting it as 0(r)( ncw ) = P(j>'(r). Fi- 
nally, we follow the path of SA from </>(r)( now ) for a length 
r. The resulting 0(0) ( now ) is the new configuration^] 



1 Note that this procedure is not inherently stable, as the one in 
[10] . but relies on the (verifiable) fact that the integration in r 
always brings sufficiently close to the saddle point. 
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The computation of the projector P is done, once for 
all, at the beginning of the simulation. However, it must 
be applied at every iteration. This can be done most 
efficiently in Fourier space, where H and P are diago- 
nal, although, for this first exploratory study on small 
lattices, we did not take advantage of this possibility. 

The cost of the algorithm depends significantly on the 
length t, that should be large enough to stretch out to 
the region where the quadratic approximation of the ac- 
tion is good. But how good is good enough? We certainly 
do not need to constrain the system on the thimble ex- 
actly, but only to the extent that the domain of integra- 
tion preserves the homology class of the thimble and the 
reweighting with the phase e lSl is feasible. 

It is then natural to ask whether t = is already 
sufficient. This corresponds to integrating the system on 
the vector space Go defined above. In general, Go does 
not belong to the same homology class as Jo, because the 
directions of steepest ascent for the quadratic part of the 
action may not, in general, be directions of convergence 
for the full action. 

However, in our simulations we observed that such di- 
vergences, although they do occur as expected, are very 
rare (see below). This suggests that the integration on 
Go, regularized, say, with a mild cutoff, might already 
provide a good approximation. Of course, such a regu- 
lator introduces an unknown bias, and the procedure is 
meaningful only if the regulator is eventually removed, by 
approaching the true thimble further. Next, we present 
our results on Go, following which we show how the true 
thimble can be systematically approached. 

Numerical results on Go — As discussed above, the 
simulations on Go are meaningful only with a regulator. 
Instead of introducing an explicit cut of the domain, we 
regularized by discarding those simulations that diverged 
within the observed histories (i.e. 4 x 10 6 trajectories for 
V = 4 4 , 10 6 trajectories for V = 6 4 and 8 x 10 5 trajecto- 
ries for V — 8 4 ). This procedure introduces an unknown 
bias, that can only be removed by approaching the thim- 
ble further. However, the fact that the divergences are 
very rare makes the regularization rather unambiguous. 
If we consider a common span of the first 8 x 10 5 trajec- 
tories, a divergence occurred with probability ~ 1.8% on 
the lattices V = 4 4 , with probability ~ 0.8% on V = 6 4 , 
and less than 0.7% on V = 8 4 (h = 5 x 10" 3 ). The 
results obtained in this way agree perfectly (within the 
rather small errors) with the results obtained with the 
algorithm of [TS] and |14j^| In particular, they show the 
correct scaling with the volume. Note that, since Go is a 
flat manifold, the residual phase, discussed in [101 111] is 
absent. 



2 We thank Gert Aarts, Christof Gattringer and Thomas Kloibor 
for sharing their partially unpublished results with us. 
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FIG. 1. Average density (n) in the critical region for the 
lattices V = 4 4 ,6 4 ,8 4 . 

We report the results of simulations for the relativistic 
Bose gas in 3 + 1 dimensions (d = 4). The mass and 
coupling were fixed at m = 1 = A, and fi was varied 
from to 1.3. In Fig. [T] and [2j we plot our results for 
the density (n) and (\(f)\ 2 ) in the most interesting range 
between \x = 0.9 and \i = 1.22. In these figures, we see a 
clear signal of the Silver Blaze phenomenon around fi ~ 
1.1. In all the simulations shown here we used At = 10~ 4 , 
but we performed also some tests with At = 10 -3 and 
At = 10 -5 and we found no significant difference. The 
errorbars on each point are computed from the standard 
deviation of 10 to 20 independent histories, in order to 
take the autocorrelation effects into account. We used 
the sources h = 5 x 10 -3 and h = 10 -3 to extract the 
limit h 0. 
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FIG. 2. Same as in Fig. [TJfor the observable (|<^| 2 ). 
In Fig. [3] we plot the average phase for the same sim- 
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FIG. 3. The data on the top-right show the average phase 
obtained with the Aurora algorithm on lattices 4 4 , 6 4 and 8 4 . 
It is interesting that the average phase is large precisely in 
the most interesting region just above /u = 1. The dashed 
lines on the bottom-left display, for comparison, the average 
phase obtained with a naive phase-quenched Monte Carlo al- 
gorithm on lattices 4 4 and 6 4 . Even on a 4 4 lattice, the sign 
problem in the phase-quenched algorithm, completely hides 
the interesting region. 



ulations reported above. The phase is used to reweight 
the observables. However, such reweighting brings cor- 
rections to the observables that are unnoticeable, within 
the statistical errors. As expected, the sign problem in Qq 
gradually increases on larger volumes and moving closer 
to the thimble will be eventually necessary. 

Moving closer to the thimble — In general, there are 
two good reasons to move closer to the thimble Jq. First, 
to remove the bias introduced by the regulator on Q . 
Second, to keep the sign problem under control on larger 
volumes. However, in the present situation, the diver- 
gences are already very rare and to observe a further 
measurable reduction would require enormous statistics. 
Moreover, the results obtained on Q are already in excel- 
lent agreement with the known results, and the reweight- 
ing with the phase has no effect even in the most critical 
case of the 8 4 lattice at [i = 1.2. Hence, the results re- 
ported here with t ^ do not intend to improve the 
precision of the results obtained above with r = 0, but 
rather to present a first exploration of the feasibility of 
moving closer to the thimble. 

To integrate Eqs. Q and ([8]), we employed the (clas- 
sical) 4th order Runge-Kutta method (RK4). This is an 
explicit method, that can be used to solve Eq. Q and Q 
as initial value problems (IVP). We argued in [10] that, 
in order to enable a stable integration in the most general 
case for large r, without the need of too tiny At, Eq. Q 
should be treated as a boundary value problem (BVP), 



by introducing explicit boundary conditions in the neigh- 
borhood of the saddle point. However, it is interesting to 
see what can be achieved even with the simpler procedure 
adopted here. 

To evaluate the closeness to the thimble, we monitored 
the reduction of the fluctuations of the imaginary part of 
the action (what really matters for the sign problem), 
when r is increased. We found that r = 4 x 10~ 2 was 
sufficient to suppress the fluctuations of the imaginary 
part of the action Si by a factor ~ 0.5 (for 4 4 ), a factor 
- 0.6 (for 6 4 ), and - 0.7 (for 8 4 ). This test was per- 
formed for fi = 1.2, in the critical region. However, a 
precise integration of the IVP becomes more and more 
difficult on increasingly large volumes (the correctness of 
the integrator can be assessed via reversibility checks). 
This shows that the IVP formulation of the ODE will 
need to be replaced by a BVP formulation in more diffi- 
cult situations. 

Finally, note that in this test we neglected the compu- 
tation of the residual phase discussed in |10|,lllj. But the 
excellent agreement with the known results, even with- 
out including the residual phase, supports the idea that 
its effect is not dramatic and maybe even negligible. 

Summary — We have reported the first numerical ap- 
plication of the Lefschetz formulation to a nontrivial 
model with a hard sign problem. In particular, we have 
studied the relativistic Bose gas model at finite chemical 
potential. Our study was restricted to small lattices, but, 
given the severity of the sign problem, this can be consid- 
ered already a very challenging test. We found excellent 
agreement with the known results already on the crud- 
est approximation of the thimble, i.e, the vector space 
Qa, once the integral was regulated by removing the few 
diverging trajectories. Moreover, we showed that it is 
possible to improve the approximation of the thimble, by 
following the equations of SD. 

Of course, the sign problem is expected to become 
worse on larger lattices: moving closer to the thimble will 
become more crucial. Work in progress include develop- 
ing efficient and stable integration algorithms to achieve 
a better approximation of the thimble, study of the scal- 
ing for larger system sizes, and application of our method 
to other models. 
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